#/usr/bin/python
# -*- coding: utf-8 -*-
'''
This module contains classes defining a domain and its components:
The class Domain represents a general one-block geometric domain while GlobalDomain 
defines the assembly of several geometric domains. Thus the class Boundary lets 
link two domains together. Class Boundary is associated with class Edge that extracts 
a geometric manifold of a domain (point, curve or surface). 

All these classes are designed in order to facilitate the implementation of domain decomposition of finite
element model by the use of primal, dual, mixte or hybrid approach [1,2]. Moreover
couples of nonconformingly meshed domains may be taken into account by the mortar methods
[3,4] through the use of Mortar Operators defined in the Boundary class. 

Details about classes are given in their own documentation string.   

References:

[1] Gosselet P., Non-overlapping domain decomposition methods in structural mechanics, 
   Rapport interne LMT-Cachan n°265, 2006.
[2] Cresta P.,  
[3] C. Bernardi, Y. Maday, and A. T. Patera, A new non conforming approach to domain 
   decomposition: The mortar element method, in Collège de France Seminar, H. Brezis
   and J.-L. Lions, eds., Pitman, 1994.
[4] Amini A., Analyse multi-échelle de structures hétérogènes par décomposition de domaine: Application aux
   navires à passagers. PhD thesis, University of Nantes, France, 2008. 


Created on 5 mars 2012

@author: Ludovic Barrière and Steven Marguet, Institut Clément Ader, 135 av. de Rangueil, 31077 Toulouse, France

© copyright:
    
    TODO
      
DocTests of the module:

'''
####### IMPORT ###########
# external library import
from mpi4py import MPI
import scipy.sparse as sp
import numpy as np
import copy as cp
# icafe library import
import sys
sys.path.append('../../../icafe')
import parameters as param

from mesh.mesh import Mesh,Constraints
# domain package library
from edge import Edge
from boundary import *
####### END IMPORT ###########


class Domain(object):
    '''
    A domain is a geometric entity that is described by one or several 
    discrete counterparts (Meshes). This class enables switching, dialoging, and so forth, from one mesh
    to another. A domain can enrich, refine or coarse its meshes.
    '''
    
    def __init__(self, userId, log=None, **kwargs):
        '''
        Domain Constructor.
        .meshesMap: Mapping between user and code meshes identifiers
        .meshes: list of meshes of the domain
        .activMesh: current mesh identifier
        .meshesLinks: Projection matrices from one mesh to another.
        
        '''
        self.UId = userId
        self.log = log
        
        self.meshes =  []
        self.meshesMap = {"U2C": {}, "C2U": []}
        self.activMesh = None
        self.edges = []
        self.edgesMap = {"U2C": {}, "C2U": []}

        # transfer operators between meshes 
        self.meshesTO = None
        self.meshesHLevel = []
        self.meshesPLevel = []
        
        # Loads and BC's:
        self.loadcases = []
        self.bC = {}
        self.lC = {}
        
        # state history:
        self.statesMap = {'STEP': [],'MODE': []}
        self.qoiMap = {'STEP': [], 'MODE': []}
        
        
        # Solution methods:
        self.analysis = None          # type of analysis
        self.incrementManager = None  # type of increment manager
        self.solvingStrategy = None   # type of solving strategy 
        
        args = self.__dict__ 
        
        for k in kwargs.keys():
            if k in args.keys():
                self.__setattr__(k,kwargs[k])
    
    def __str__(self):
        
        r = '\n' + 'Domain ' + str(self.UId) + ':\n'
        r += '  - meshes: ' + str(self.meshesMap['C2U']) + ' (mesh ' + str(self.activMesh) + ' activated)'+'\n'
        r += '  - edges:  ' + str(self.edgesMap['C2U']) + '\n'
        r += '  - current step: ' + str(self.currentStep)
        return r  
    
    def getUId(self):
        
        return self.UId
     
    ## LOAD AND BOUNDARY CONDITIONS:  
    def initStep(self, stepId=None):
        '''
            Initialize a step in domain
            @param stepId: str, optional. step identifier, for information only.
        '''
        self.initLoadCases()
        self.currentStep = stepId
        
        for mesh in self.meshes:
            mesh.initConstraints()
        
        self.saveFields('ST'+str(stepId)+'INIT', level='step')
        
        return 0
          
    def initLoadCases(self):
        '''
            Initializes load and boundary conditions and load cases.
        '''
        self.loadcases = []
        self.bC = {}
        self.lC = {}
            
    def addLoadCase(self, loadcase):
        '''
            Adds a load case.
        '''
        self.loadcases.append(loadcase)
        
    def rmLoadCase(self, loadcase):
        '''
            Removes a load case.
        '''
        self.loadcases.remove(loadcase)
        
    def addBC(self, BC):
        '''
            Adds a boundary condition.
        '''
        self.bC[BC.UId] = BC
    
    def rmBC(self, BC):
        '''
            Removes a boundary condition.
        '''
        try:
            del(self.bC[BC.UId])
        except:
            self.log.warning("The boundary condition to remove was not found in the domain" + str(self.UId))    
        
    def addLoad(self, load):
        '''
            Adds a load condition.
        '''
        self.lC[load.UId]=load
    
    def rmLoad(self, load):
        '''
            Removes a load condition.
        '''
        try:
            del(self.lC[load.UId])
        except:
            self.log.warning("The load condition to remove was not found in the domain" + str(self.UId))
     
        
    ## OPERATORS: 
    
    def setConstraints(self, meshId=None, stepId=None):
        '''
            Sets and returns the constraints on primal (displacements) and dual (forces) degrees of freedom.
            Usage:
            Uprescr[i] = Cd[i].T * d[i]
            Fext[i] = Cf[i].T * f[i]
            i being the load case number.
            
            @param meshId: int, optional. 
            @param stepId: str, optional.
            
            @return Cd: list of sparse matrix. Extractors of prescribed displacements (one per load case).
            @return Cf: list of sparse matrix. Extractor of external forces (one per load case).
            @return d: list of sparse matrix. Prescribed values of displacements (one per load case).
            @return f: list of sparse matrix. Values of external forces (one per load case). 
        ''' 
        if not meshId:
            meshId=self.activMesh
        
        # Get mesh info:
        mesh = self.meshes[meshId]
        
        return mesh.setConstraints(self.loadcases, self.bC, self.lC)
        
        
        
        
        mesh.setConstraints(constraints)
            
        return constraints
    
    def getConstraints(self, meshId=None, full=True):
        '''
            Direct access to constraints built by the setConstraints() method. If not existing in mesh, calls the setConstraints() method.
            @param meshId: str or int. Identifier of the mesh. Optional.
            
            @return Cd
            @return d
            @return Cf
            @return f
        '''
        if not meshId:
            meshId=self.activMesh
        
        constraints = self.meshes[meshId].getConstraints()
        
        if not constraints:
            constraints = self.setConstraints(meshId)
        
        return constraints
        
    def getLhs(self, lhs, meshId=None):
        if not meshId:
            meshId=self.activMesh
        if isinstance(lhs,list):
            lhs = [self.meshes[meshId].getLhs(lhs_name) for lhs_name in lhs]
        else:
            lhs = self.meshes[meshId].getLhs(lhs)
        
        return lhs
        
    def getRhs(self, rhs, meshId=None):
        if not meshId:
            meshId=self.activMesh
        if isinstance(rhs,list):
            rhs = [self.meshes[meshId].getRhs(rhs_name) for rhs_name in rhs]
        else:
            rhs = self.meshes[meshId].getRhs(rhs)
        
        return rhs
    
    def setLhs(self, lhs, lhs_name):
        ''' 
            Not needed actually because operators are stored in meshes instances.
        '''
        return NotImplementedError
    
    def setRhs(self, rhs, rhs_name):
        ''' 
            Not needed actually because operators are stored in meshes instances.
        '''
        return NotImplementedError 
    
    
    def computeOperators(self, lhslist, rhslist, disp=1, load=1, meshId=None):
        
        if not meshId:
            meshId=self.activMesh
            
        lhs, rhs = self.meshes[meshId].computeOperators(lhslist, rhslist, disp, load)
        
        return lhs, rhs
    ## MULTI-MESH MANAGEMENT:
    
    def upDateDofs(self, ddofs, dduals, meshId=None, op='add'):
        
        if not meshId:
            meshId=self.activMesh
            
        self.meshes[meshId].upDateDofs(ddofs, dduals, op)
        
        
        return self.meshes[meshId].getDofs(format=False)
                
    
    REMESH_MODES = {
        "default": "Replace",
        "R": "Replace", 
        "AG": "Add Grid",
        }
    
    def hRefineMesh(self, MODE = "default"):
        return NotImplementedError
        
    def pRefineMesh(self, MODE = "default"):
        return NotImplementedError
    
    def hpAdaptMesh(self, MODE = "default"):
        return NotImplementedError
    
    def addMesh(self, mesh):
        print 'adding mesh in domain ', self.UId
        codeId = len(self.meshesMap["C2U"])
        self.meshes.append(mesh)
        if mesh.UId not in self.meshesMap["C2U"]:
            self.meshesMap["U2C"][mesh.UId] = codeId
            self.meshesMap["C2U"].append(mesh.UId) 
        else:
            return KeyError("Mesh " + repr(mesh.UId) + " already exists.")
        if codeId == 0:
            self.activMesh = 0
        
        # add mesh to edges:
        for edge in self.edges:
            edge.addMesh(mesh)
        
        return codeId
    
    def setActivMesh(self, meshCId = None, meshUId = None):
        '''
            Sets the active mesh of the domain.
            @param meshCId: int, optional.
            @param meshUId: int, optional.
        '''
        
        if meshCId:
            self.activMesh = meshCId
        if meshUId:
            self.activMesh = self.meshesMap["U2C"][meshUId]
    
    def getActivMesh(self):
        '''
            Returns the active mesh of the domain.
        '''
        return self.meshes[self.activMesh]
    
    
    def buildTransferOperators(self, mesh1CId, mesh2CId):
        '''
            Computes the transfer operators between two meshes.
            @param mesh1CId: int.
            @param mesh2CId: int.
        '''
        return NotImplementedError
    
    
    ## OTHER STUFFS:
    
    def addEdge(self, UId, edge=None, nodelistCId=None, nodelistUId=None, meshUId=None, point_list=None, coord_expression=None, error=10e-6):          
        '''
          Adds edge to the Domain.
          @param edge is an Edge object to add to the domain.
          @param nodelist is a list of nodes defining the edge.
          @param meshUId is the user id of the mesh from which the nodes come. 
          @param point_list is a list of point coordinates that define an edge.
          @param coord_expression is a tuple containing an analytical expression of coordinate relation. 
          For instance, the following coord_expression "[{x},2*{x},0]" matches points on the line y=2x and z=0.
          @param error defines the approximation error to the coordinate of points. default value is 10e-6.
        '''
        try:
            meshCId=self.meshesMap["U2C"][meshUId]
        except KeyError:
            meshCId = self.activMesh

        if edge:
            # edge is given
            pass
        elif nodelistCId or nodelistUId:
            # build the edge from a node list and a mesh.
            edge = Edge( UId, self.UId, nodelistCId=nodelistCId, nodelistUId=nodelistUId, mesh=self.meshes[meshCId])
        elif point_list:
            # build the edge from a point list.
            edge = Edge(UId, self.UId,pointlist=pointlist, mesh=self.meshes[meshCId])
        elif coord_expression:
            # build the edge from a coordinate equation.    
            edge = Edge(UId, self.UId, coord_expression=coord_expression, mesh=self.meshes[meshCId])
            
        self.edges.append(edge)    
        self.edgesMap["U2C"][edge.UId]=len(self.edges)-1
        self.edgesMap["C2U"].append(edge.UId)
        
        return edge
    
    def getNodesFromEdge(self, edge, mesh=None):
        '''
            Return a list of nodes belonging to the edge
            @param edge: Instance of class Edge.
        '''
        # init list of nodes:
        nodes = []
        
        if not mesh:
            mesh = self.meshes[self.activMesh]
        
        try:
            eid = self.edgesMap['U2C'][edge.UId]
            nodes = self.edges[eid].getNodesFromMesh(mesh)
            
        except KeyError: # if edge does not correspond to an edge of the domain:
            nodes = []
        
        return nodes
    
    def saveResults(self):
        '''
            A function to save the resulting fields of the analysis. 
        '''
        pass
    
    def saveFields(self, stateKey, fields=None, level='iter', meshes=None, type=None):
        '''
           A function to save the fields of the current state.
           @param fieldKey: keyword of the state
           @param fields: tuple of fields identifiers to be saved, default=('dof', 'stress', 'strain',)
           @param meshes: list of meshes which will save the retained fields, default=None (i.e. domain.activMesh)
           @param type: save fields as particular state in a equilibrium time/load step (type=STEP) or as a unreferenced mode (type=MODE) 
        '''
        if not stateKey:
            return 0
        
        if not meshes: meshes=[self.activMesh]
        
        if not type:
            if level=='mode':
                type = 'MODE'
            else:
                type = 'STEP'
                
        # check the key:
        try:
            self.statesMap[type].index(stateKey)
            #raise AttributeError('Keyword of the field has already been addressed: ', stateKey)
        except:
            self.statesMap[type].append(stateKey)
        
        if not fields:
            param = parameters.load()
            fields, components = param.getOutputFields(level)
        
        # save the field in meshes:
        for m in meshes:
            self.meshes[m].saveFields(stateKey, fields, components)
            
    def saveQOI(self, stateKey, qoi=None, level='iter', meshes=None, type=None):
        
        if not stateKey:
            return 0
        
        if not meshes: meshes=[self.activMesh]
        
        if not type:
            if level=='mode':
                type = 'MODE'
            else:
                type = 'STEP'
                
        try:
            self.qoiMap[type].index(stateKey)
        except:
            self.qoiMap[type].append(stateKey)
            
        if not qoi:
            param = parameters.load()
            qoi, args = param.getOutputQOI(level)
        else:
            args = tuple([None]*len(qoi)) 
                 
        for m in meshes:
            self.meshes[m].saveQOI(stateKey, qoi, args)
                        
    def setAnalysis(self, analysis):
        '''
           sets the instance of analisys.Analysis() as attribute of the domain
        '''
        self.analysis = analysis
        
    def setIncrementMngr(self, incrementMngr):
        '''
        sets the instance of stepManager.IncrementManager() as attribute of the domain
        '''
        self.incrementManager = incrementMngr 
        
    def setSolvingStrategy(self, strategy):
        '''
        sets the instance of strategy.SolvingStrategy() as attribute of the domain
        '''
        self.solvingStrategy = strategy
    
    def analyse(self, meshId=None, updateOperators=True):
        '''
            Is the self/analysing method of domain which uses instances of Analysis(), IncrementManager(), Solvingstrategy(), etc.
        '''
        
        lhslist, rhslist = self.analysis.getOperators()
        
        disp, load = self.incrementManager.setNextIncr()
        
        print 'initial increment : disp=', disp, ' load=', load
        
        while 1: # increment loop
            
            if self.analysis.iterative:
                print 'solve a non-linear problem'
                # solve nonlinear problem:
                while 1: # iteration loop
                    
                    # compute operators and residual    
                    if self.incrementManager.updateLhs:
                        cur_lhslist = copy(lhslist)
                    else:
                        cur_lhslist = []                
                        
                    res = self.solvingStrategy.computeOperators(self, cur_lhslist, rhslist, disp, load, meshId, error=True)
                    
                    # convergence criterion:
                    if res < self.analysis.tol: break 
                    
                    ddofs, dduals = self.solvingStrategy.solve(self, lhslist, rhslist, disp, load, meshId)
                    
                    # prepare next iteration and compute continuation parameters:
                    stateKey = self.incrementManager.getStateKey(level='iter')
                    load, ddofs = self.incrementManager.setNextIter(ddofs)
                     
                    # increment d.o.f.:
                    self.upDateDofs(ddofs, dduals, meshId)
                    
                    self.saveFields(stateKey, level='iter')
                    self.saveQOI(stateKey, level='iter')
                    
                    if not disp and not load: 
                        break # end of iterations for poor convergence reasons
                    
                    
#            elif self.analysis.modal:
#                print 'solve a modal problem'
#                # compute operators and residual                    
#                self.solvingStrategy.computeOperators(self, lhslist, rhslist, disp, load, meshId, error=False)
#                
#                # solve linear problem:
#                ddofs = self.solvingStrategy.solve(self, lhslist, rhslist, meshId)
#                
#                # update dof:
#                self.upDateDofs(ddofs, meshId)
#            
#                stateKey = self.incrementManager.getStateKey(level='mode')
#                self.saveFields(stateKey)
#                self.saveQOI(stateKey)
#                
            else:
                print 'solve a linear problem'
                print 'disp = ', disp, ' load = ', load
                # compute operators and residual                    
                self.solvingStrategy.computeOperators(self, lhslist, rhslist, disp, load, meshId, error=False)
                
                # solve linear problem:
                ddofs, dduals = self.solvingStrategy.solve(self, lhslist, rhslist, meshId)
                
                # update dof:
                self.upDateDofs(ddofs, dduals, meshId, op='replace')
                
                stateKey = self.incrementManager.getStateKey(level='incr')
                self.saveFields(stateKey, level='incr')
                self.saveQOI(stateKey, level='incr')

#                print ddofs.todense()
                
            # a posteriori error measurment
            # ...
            # save increment solution:
            
            # set next increment:
            disp, load = self.incrementManager.setNextIncr()
            if not disp and not load: break # end of increments
        
        return 0  
        
        
        
class DomainDecorator(Domain):                
    ''' Minimal dynamic decorator of class Domain.'''
    
    def __init__(self,domain, intracomm=None, parent_intercomm=None, child_intercomm=None):
        '''Constructor of the decorator. '''
        kwargs=domain.__dict__
        super(DomainDecorator,self).__init__(kwargs["UId"],**kwargs)
        #self.domain = domain
        self.master = parent_intercomm
        self.coworkers = intracomm
        self.workers = child_intercomm
        self.jobs = {0: self.recv,
                1: self.send,
                2: self.reduce, 
                3: self.scatter,
                5: self.isend,
                6: self.irecv,
        }
        
    def waitForJob(self, parent_intercomm=None, intracomm=None):
        
        msg=None
        
        while 1:
            jobid,args,kwargs = self.master.bcast(root=0)
            
            if jobid == 'staticmethod':
                args[0](self, *args[1:], **kwargs)
            
            elif jobid=='disconnect':
                try:
                    self.workers.bcast(('disconnect',[],{}),root=MPI.ROOT)
                except AttributeError:
                    pass
                msg='disconnect'
                break
            
            elif jobid=='domain':
                msg='domain'    
            
            try:
                jobid = self.jobs[jobid](*args,**kwargs)
            except KeyError:
                msg = self.call_method(jobid,*args,**kwargs)
                    
            if msg:
               break 
            
        return msg
    
    def call_method(self, method, *args, **kwargs):
        
        try:
            m = getattr(self,method)
        except AttributeError:
            pass
        else:
            return m(*args,**kwargs)            
                 
    def recv(self, parent_intercomm=None):   
        pass
    
    def send(self, parent_intercomm=None):
        pass
    
    def reduce(self):
        pass
    
    def scatter(self):
        pass
    
    def irecv(self):
        pass
    
    def isend(self):
        pass
    
    
    
     
class GlobalDomain(DomainDecorator):
    '''
    The class GlobalDomain() inherits from the class Domain(). It implements a decorator for 
    instances of class Domain(). 
    '''
    
    def __init__(self, domain, subDomains=[]):
        '''
        GlobalDomain Constructor:
        1. Constitutive domains definition:
        .domains: list of constitutive domains
        .domainsMap: Mapping between user and code identifiers of the domains
        .domainsEdges: lists of edges of each domain
        .domainsLhs: dictionary of Left-Hand-Sides of equilibrium equation(s)
        .domainsRhs: dictionary of Right-Hand-Sides of equilibrium equation(s)
        
        2. Boundaries definition:
        .boundaries: list of domains and corresponding edges forming boundaries between two domains
        .boundariesMap: Mapping between user and code identifiers of the boundaries
        .boundariesTO: list of trace operators from associated domains towards the boundary
        .boundariesBAO: list of binary assembly operators from domains towards global domain 
        .boundariesSAO: list of signed assembly operators from domains towards global domain
        .boundariesType: list of type of boundaries ( 0 if primal, 1 if dual and 2 if mixed)
        .boundariesConstraints: list of constraints to be applied on d.o.f. of boundaries ( C q = d )
        .boundariesScaling: list of scaling matrices
        
        '''
        # initialize domain:
        super(GlobalDomain,self).__init__(domain)
        # initialize subdomains:
        self.domains = []
        self.domainsMap = {"U2C": {}, "C2U": []}
#        self.domainsEdges = []
        self.domainsLhs = {}
        self.domainsRhs = {}
        self.domainsBoxes = []
        self.domainsConnect = {}
        
        self.boundaries = []
        self.boundariesMap = {"U2C": {}, "C2U": []}
        self.boundariesTO = []
        self.boundariesBAO = []
        self.boundariesSAO = []
        self.boundariesTypes = []
        self.boundariesGI = []
        self.boundariesInterpMesh = []
        
        self.bdofs = 0
#        self.boundariesConstraints = []
#        self.boundariesScaling = []
        
        if len(subDomains) > 0:
            for i, d in enumerate(subDomains):
                self.addDomain(d)  
        else:
            # domains are predefined in node sets and element sets:
            pass
        
        if len(subDomains)>1:
            self.buildLocalBoundaries()
        
        if param.PARALLEL:
            self.spawnProcesses(data='full')
        else:
            self.workers=None
            
    def spawnProcesses(self, data='light'):
        # create processes:
        self.workers = MPI.COMM_SELF.Spawn(sys.executable,['workers.py'],maxprocs=len(self.domains))
        
        # initialize subdomains:
        if data=='light':
            for domain in self.domains:
                domain.purgeData()
        
        requests = []
        for i,di in enumerate(self.domains):
            requests.append(self.workers.isend(di,dest=i))
        MPI.Request.Waitall(requests)
#        for domain in self.domains:
#            self.workers.isend(di,dest=i,root=MPI.ROOT)
        self.domains = [None]*len(self.domains)
           
    def unSpawnProcesses(self, data='light'):
        
        if data:
            # retrieve subdomains:
            self.workers.bcast(('domain',[],{}),root=MPI.ROOT)
            self.domains = self.workers.gather(self.domains,root=MPI.ROOT)
        
        self.workers.bcast(('disconnect',[],{}),root=MPI.ROOT)
        
        self.workers.Disconnect()
       
    # domain methods overload:
    def getLhs(self, lhs, meshId=None, full=False):
        
        if isinstance(lhs,list):
            out = [None]*len(lhs)
            L = True
        else:
            out = None
            L = False
        
        for domain in self.domains:    
            lhs_d = domain.getLhs(lhs,meshId)
            if L:
                try:
                    out = [out[i] + lhs_d[i] for i in len(out)]
                except:
                    out = lhs_d
            else:
                try:
                    out += lhs_d
                except:
                    out = lhs_d
        
        return out
    
    def getRhs(self, rhs, meshId=None, full=False):
        
        if isinstance(rhs,list):
            out = [None]*len(rhs)
            L = True
        else:
            out = None
            L = False
        
        for domain in self.domains:    
            rhs_d = domain.getRhs(rhs,meshId)
            if L:
                try:
                    out = [out[i] + rhs_d[i] for i in len(out)]
                except:
                    out = rhs_d
            else:
                try:
                    out += rhs_d
                except:
                    out = rhs_d
        
        return out
    
    def setConstraints(self, meshId=None, stepId=None):
        '''
            Sets the constraints in subdomains.
        '''
        
        
    
    def getConstraints(self, meshId=None, full=False):
        '''
            Gets the constraints on boundary or on the entire subdomain.
            @param meshId:
            @param full: True or False, optional. 
        '''
        
        
        
        return constraints
        
    def computeOperators(self, lhslist, rhslist, updateDofs=False, updateInputFile=True, runThirdSoftware=True, readOperators=True):
        
        if param.PARALLEL:
            args = [lhslist,rhslist]
            kwargs = {'updateDofs': updateDofs, 'updateInputFile': updateInputFile, 'runThirdSoftware': runThirdSoftware,
                       'readOperators': readOperators}
            self.workers.bcast(('computeOperators',args,kwargs),root=MPI.ROOT)
        else:
            for domain in self.domains:
                domain.computeOperators(lhslist,rhslist,updateDofs)
        
    # own methods:    
    def splitDomain(self):
        pass
        
    def joinDomains(self):
        pass
        
    def addDomain(self, domain):
        
        codeId = len(self.domainsMap["C2U"])
        if  not str(domain.UId).startswith(str(self.UId)+'.'):
            domain.UId = str(self.UId)+'.'+str(domain.UId) 
        self.domains.append(domain)
        if domain.UId not in self.domainsMap["C2U"]:
            self.domainsMap["U2C"][domain.UId] = codeId
            self.domainsMap["C2U"].append(domain.UId)
            #self.domainsEdges.append([])          
        else:
            return KeyError("Domain " + repr(domain.UId) + " already exists.")
    
    def buildLocalBoundaries(self):
        '''
           Creates local boundaries (i.e. boundaries between two domains).
           The method consists in finding shared nodes and building edges for each domain with them. 
           The two edges are then associated in a LocalBoundary instance.
        '''
        
        # find neighbouring subdomains.
        self.buildBoxes()
        self.buildConnectivityOfDomains()
        
        for did in xrange(len(self.domains)):
            for neighbId,sharedItems in self.domainsConnect[did].iteritems():
                
                if neighbId > did:
                # find shared nodes of connected subdomains.
                    if not sharedItems:
                        sharedItems = self.findSharedItems(did,neighbId)
                    # build edges
                    if sharedItems:
                        edge1 = self.addGlobalEdge(did, sharedItems = sharedItems)
                        try:
                            edge2 = self.addGlobalEdge(neighbId, sharedItems = self.domainsConnect[neighbId][did])
                        except KeyError:
                            raise KeyError("Domains' connectivity is corrupted")
                        # associated edges in a local boundary
                        print 'add boundary between domains ', did, ' and ', neighbId
                        self.addLocalBoundary(did, len(self.edges)-2, neighbId, len(self.edges)-1)
                        
                        assert self.domains[did].meshes[0].nodesMap['C2COO'][self.domainsConnect[did][neighbId][0]]==self.domains[neighbId].meshes[0].nodesMap['C2COO'][self.domainsConnect[neighbId][did][0]]
#                        print 'node 0 of boundary:', self.domains[did].meshes[0].nodesMap['C2COO'][self.domainsConnect[did][neighbId][0]]
#                        print 'node 1 of boundary:', self.domains[did].meshes[0].nodesMap['C2COO'][self.domainsConnect[did][neighbId][1]]
#                        print 'node 2 of boundary:', self.domains[did].meshes[0].nodesMap['C2COO'][self.domainsConnect[did][neighbId][2]]
                    else:
                        # not connected domains: remove from connectivity table.
                        self.domainsConnect[did].pop(neighbId)
                        self.domainsConnect[neighbId].pop(did)
        
    
    def buildConnectivityOfDomains(self):
        '''
            Build the domainsConnect dictionary with . Key of the dictionary is domain's code identifier. Value is a dictionary of 
            neighboring domain's id as keys and shared items as values:
            For example, 3 domains (0 connected with 1 and 2, 1 and 2 are not connected with each other) with empty shared items (None):
            domainsConnect = { 0: {1:None, 2:None}, 1: {0:None}, 2:{0:None} }
        '''
        self.domainsConnect = {}
        
        for did1 in xrange(len(self.domains)):
            for did2 in xrange(did1+1,len(self.domains)): 
                if not self._are_disconnected(did1,did2):
                    #print did1, did2
                    self.domainsConnect[did1]={}
                    self.domainsConnect[did1][did2] = None
                    self.domainsConnect[did2]={}
                    self.domainsConnect[did2][did1] = None
    
    def buildBoxes(self):    
        '''
            build the 3D box which fits domains :
            box = np.array([[xmin, xmax],
                           [ymin, ymax], 
                           [zmin, zmax]])
        '''
        self.domainsBoxes = []
        for d in self.domains:
            box = np.zeros((3,2),dtype = 'd')
            nodescoo = d.meshes[0].getNodesArray(output_type=np.array)
            box[0,0] = np.min(nodescoo[:,0])
            box[0,1] = np.max(nodescoo[:,0])
            box[1,0] = np.min(nodescoo[:,1])
            box[1,1] = np.max(nodescoo[:,1])
            box[2,0] = np.min(nodescoo[:,2])
            box[2,1] = np.max(nodescoo[:,2])
            self.domainsBoxes.append(box)
            
            
    def _are_disconnected(self, did1, did2):
        
        disconnected = False
        
        try:
            box1 = self.domainsBoxes[did1]
            box2 = self.domainsBoxes[did2]
        except IndexError:
            self.buildBoxes()
            box1 = self.domainsBoxes[did1]
            box2 = self.domainsBoxes[did2]
            
            
        for i in xrange(3):
            xmin = max([box1[i,0],box2[i,0]])
            xmax = min([box1[i,1],box2[i,1]])
            if xmax < xmin:
                disconnected = True
                break
        
        return disconnected
         
    def findSharedItems(self,did1, did2, method='matching'):
        '''
        
        @param method: str, optional. Defines the method used to look for shared items. Default value is 'matching' that means
        shared items are looked for in matching nodes only. 'non-matching' enables to get nodes of one domain that does not 
        match any node of second domain whereas they lie on its interpolated geometry.
        '''
        sharedItems = [None,None]
        
        if method == 'matching':
            for nid1,node in enumerate(self.domains[did1].meshes[0].nodesMap["C2COO"]):
                try:
                    nid2 = self.domains[did2].meshes[0].nodesMap["COO2C"][node]
                except KeyError:
                    continue
                
                if not sharedItems[0]:
                    sharedItems[0]=[]
                if not sharedItems[1]:
                    sharedItems[1]=[]
                    
                sharedItems[0].append(nid1)
                sharedItems[1].append(nid2)
                
                
        if not None in sharedItems:
            self.domainsConnect[did1][did2]=sharedItems[0]
            self.domainsConnect[did2][did1]=sharedItems[1]
        
        return sharedItems[0]
            
    def addLocalBoundary(self, domain1CId, edge1CId, domain2CId=None, edge2CId=None, boundaryUId = None):
        '''
        Adds a local boundary defined by one domain at least (two domains at most) and corresponding edge(s).
        @param domain1CId: integer, identifier of the first domain of the local boundary.
        @param edge1CId: integer, code identifier of the edge of the first domain.
        @param domain2CId: integer, optional. Code identifier of the second domain of the local boundary.
        @param edge2CId: integer, optional. Code identifier of the edge of the second domain.
        @param boundaryUId: str or int, optional. User identifier of the local boundary.
        '''
        if not boundaryUId:
            if domain2CId:
                boundaryUId = 'B.' + str(self.edges[edge1CId].UId)+'-'+str(self.edges[edge2CId].UId)
            else:
                boundaryUId = 'B.' + str(self.edges[edge1CId].UId)+'-outside'
            
        nboundary = len(self.boundaries)
        self.boundaries.append([domain1CId,edge1CId,domain2CId,edge2CId])
        self.boundariesMap["C2U"].append(boundaryUId)
        self.boundariesMap["U2C"][boundaryUId]=nboundary
        self.boundariesTypes.append('P')
        
        mesh1UId = self.domains[domain1CId].getActivMesh().UId
        if domain2CId:
            mesh2UId = self.domains[domain2CId].getActivMesh().UId
        
        edge1mesh = self.edges[edge1CId].getInterpMesh(mesh1UId)
        if domain2CId:
            edge2mesh = self.edges[edge2CId].getInterpMesh(mesh2UId)
        else:
            edge2mesh=None
            
        boundarymesh, mortars = LocalBoundary.setInterpMesh(edge1mesh,edge2mesh)
        
        self.boundariesInterpMesh.append(boundarymesh)
                                         
        BAO, SAO = LocalBoundary.buildAssemblyOperators(boundarymesh,edge1mesh,edge2mesh,mortars)
        self.boundariesBAO.append(BAO)
        self.boundariesSAO.append(SAO)
        self.boundariesTO.append([self.edges[edge1CId].getTraceOperator(mesh1UId), self.edges[edge2CId].getTraceOperator(mesh2UId)])
        self.boundariesGI[i][0]=self.bdofs
        if self.boundariesTypes[i] == 'M':
            self.bdofs += self.boundariesBAO[i][0].get_shape()[0] 
        self.boundariesGI[i][1]=self.bdofs
        self.bdofs += self.boundariesBAO[i][0].get_shape()[0]
    
    def ResetGlobalIndex(self):
        '''
            Computes the index of the first d.o.f. of the local boundaries in the global boundary d.o.f. numbering. 
        '''
        # compute the total number of boundary d.o.f.:
        self.bdofs = 0
        for i,bd in enumerate(self.boundaries):
            self.boundariesGI[i][0]=self.bdofs
            if self.boundariesTypes[i] == 'M':
                self.bdofs += self.boundariesBAO[i][0].get_shape()[0] 
            self.boundariesGI[i][1]=self.bdofs
            self.bdofs += self.boundariesBAO[i][0].get_shape()[0] 
            
            
#        # compute the global assembly operators of each boundary:
#        bdofs = 0
#        self.boundariesBAOg = []
#        self.boundariesSAOg = []
#        for i,bd in enumerate(self.boundaries):
#            dim1 = self.boundariesBAO[i][0].get_shape()
#            dim2 = self.boundariesBAO[i][1].get_shape()
#            
#            assert dim1[0]==dim2[0], "Shape of local assembly operators mismatch"
#            
#            BAOg = [sp.lil_matrix((self.bdofs,dim1[1])),sp.lil_matrix((self.bdofs,dim2[1]))]
#            SAOg = [sp.lil_matrix((self.bdofs,dim1[1])),sp.lil_matrix((self.bdofs,dim2[1]))]
#            
#            # domain 1:
#            BAOg[0][bdofs:bdofs+dim1[0],:] = self.boundariesBAO[i][0]
#            SAOg[0][bdofs:bdofs+dim1[0],:] = self.boundariesSAO[i][0]
#            
#            # domain 2:
#            BAOg[1][bdofs:bdofs+dim2[0],:] = self.boundariesBAO[i][1]
#            SAOg[1][bdofs:bdofs+dim2[0],:] = self.boundariesSAO[i][1]
#            
#            bdofs += dim2[0]
#            
#            self.boundariesBAOg.append(BAOg)
#            self.boundariesSAOg.append(SAOg)  
#        
#        assert bdofs == self.bdofs, "Global assembly operators are not valid"
        
    def initSubDomainsConnection(self):
        '''
            Initializes the connections of subdomains by passing assembly and trace operators down to them, for each boundary to which 
            they belong.            
        '''
        
        
         
        
    def addGlobalEdge(self, domainId, sharedItems):          
        '''
          Adds edges to the GlobalDomain depending on the number of coincident subdomains.
          @param point_list is a list of point coordinates that define an edge.
          @param coord_expression is a tuple containing an analytical expression of coordinate relation. 
          For instance, the following coord_expression "[{x},2*{x},0]" matches points on the line y=2x and z=0.
          @param error defines the approximation error to the coordinate of points. default value is 10e-6.
        '''
        UId = str(domainId) + '.'
        # make the edge user id unique:
        i=0
        while 1:
            if UId + str(i) in self.domains[domainId].edgesMap["C2U"]:
                i+=1
            else:
                UId += str(i)
                break
               
        newedge = self.domains[domainId].addEdge(UId, nodelistCId=sharedItems)
        
        self.addEdge(newedge.UId, edge=newedge)
        
        return newedge
        
class SubDomain(DomainDecorator):
    
    def __init__(self, domain, intracomm=None, intercomm=None):
        
        # initialize domain:
        super(SubDomain,self).__init__(domain)
        
        # define a container of boundaries operators and types: 
        self.connections = {} # keys: UId of connected subdomain, values: list of type, trace and assembly operators of the boundary
        
        
#    Unchanged functions from the class Domain:        
#    def setConstraints(self, meshId=None, stepId=None):
#        '''
#            Sets the constraints on the subdomain.
#        '''
#        pass
#    
#    def getConstraints(self, meshId=None):
#        '''
#            Gets the constraints on the subdomain.
#        '''
#        pass    
    
    def setTraceOperators(self, traceOperators):
        self.traceOperators = traceOperators
    
    def setAssemblyOperators(self):
        pass
        
    def addConnection(self,domainUId,type,Top,Aop):
        self.connections[domainUId] = [type,Top,Aop]
        
    def getCondensedLhs(self, lhslist, traceOperators=None, meshId=None, output_matrix=sp.lil_matrix):
        '''
            Returns the condensed lhs operators of the subdomain's mesh.
            @param traceOperators: tuple.
            @param lhslist: list of lhs, optional. Default value is None.
            @param meshId: Identifier of the mesh of the subdomain, optional. Default value is None.
            @param output_matrix: sparse matrix class, optional. Specifies the parent class of the output.  
            @return S: sparse matrix in the output_matrix format. Hybrid Schur complement built from the Lhs.
        '''
        
        if traceOperators:
            self.setTraceOperators(traceOperators)
        elif not traceOperators and not self.traceOperators:
            raise ValueError("No trace operators were given for the building of condensed LHS operators of subdomain " + str(self.UId))
        
        TOp = []
        TOd = []
        TOm = []
        TOi = None
        # treatment of trace operators:
        try:    
            TOp, TOd, TOm = self.traceOperators
        except ValueError:
            if len(self.traceOperators==4):
                TOp, TOd, TOm, TOi = self.traceOperators
            else:
                raise ValueError("Trace Operators of subdomain " + str(self.UId) + " is not correct.")
        
        TO=[None, None, None]
        for t in TOp:
            try:
                TO[0] = sp.bmat([[TO[0]],[t]])
            except:
                TO[0] = t
        for t in TOd:
            try:
                TO[1] = sp.bmat([[TO[1]],[t]])
            except:
                TO[1] = t
        for t in TOm:       
            try:
                TO[2] = sp.bmat([[TO[2]],[t]])
            except:
                TO[2] = t
        
        # build trace operator of inner d.o.f.
#        if not TOi:
#            TOi 
        
        P = True
        D = False
        M = False
        
        
        ## Init operators:
        # trace operators:
        TO = [TO[0].todense(), TO[1].todense(), TO[2].todense()]

        TOi = TOi.todense()
        TOx = sp.bmat([[TOi],[TO[1]],[TO[2]]]).todense()
        TOd
        TOm
        
        
        # assembly operators:
        Ap       # Primal assembly operator
        Ad       # Dual assembly operator
        
        # Not Implemented: Mixed boundary operators
        #Am = None  # Mixed assembly operator
        #Am_n       # Mixed assembly operator of the neighbouring subdomains
        #Tm         # Boundary equivalent stiffness of the subdomains
        #Tm_n       # Boundary equivalent stiffness from the neighbouring subdomains
        
        S = []
        
        for lhs in self.getLhs(lhslist, meshId):
            lhs = lhs.todense()
            
            lhs_xx = TOx.dot(lhs).dot(TOx.T)
            lhs_xx_inv = npla.pinv(lhs_xx)
            lhs_px = TOp.dot(lhs).dot(TOi.T)
            lhs_pp = TOp.dot(lhs).dot(TOp.T)
            
            # primal part
            if P:
                Sp = lhs_pp - lhs_px.dot(lhs_xx_inv).dot(lhs_px.T)
                
            # dual part
            if D:
                Sd = TOd.dot(TOx.T).dot(lhs_xx_inv).dot(TOx).dot(TOd.T)
            
            # mixed part
            if M:
                Sm = None # TOm.dot(TOx.T).dot(lhs_xx_inv).dot(TOx).dot(TOm.T)
            
            # coupling parts
            if P and D:
                Spd = lhs_px.dot(lhs_xx_inv).dot(TOx).dot(TOd.T)
            if P and M:
                Spm = None
            if D and M:
                Sdm = None
            
            # assembling parts:
            if P and D and M:
                Spdm = sp.bmat([[Sp,Spd,Spm],[Spd.T,Sd,Sdm],[Spm.T,Sdm.T,Sm]])
            elif P and D:
                Spdm = sp.bmat([[Sp,Spd],[Spd.T,Sd]])
            elif P and M:
                Spdm = sp.bmat([[Sp,Spm],[Spm.T,Sm]])
            elif D and M:
                Spdm = sp.bmat([[Sd,Sdm],[Sdm.T,Sm]])
            
            S.append(Spdm)
            
        return S

    def getCondensedRhs(self, rhslist, traceOperators=None, ):
        pass
    
    def localization(self):   
        pass
        

        
        
        
           
        